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ABSTRACT 

We  examine  the  problem  of  maximum  likelihood  covari¬ 
ance  estimation  using  a  sensor  array  in  which  the  relative 
positions  of  individual  sensors  change  over  the  observa¬ 
tion  interval  The  problem  is  cast  as  one  of  estimating 
a  structured  covariance  matrix  sequence.  A  vector  space 
structure  is  imposed  on  such  sequences,  and  within  that 
vector  space  we  define  a  constraint  space  given  by  the  in¬ 
tersection  of  a  hyperplane  and  the  space  of  sequences 
of  nonnegative  definite  matrices  W2.  Knowledge  of  the 
changing  array  geometry  is  used  to  reduce  the  dimension 
of  the  search  space.  An  extension  of  the  inverse  iteration 
algorithm  of  Burg  etal  is  proposed  for  finding  the  maxi¬ 
mum  likelihood  solution. 

1.  INTRODUCTION 

In  many  array  signal  processing  applications  knowledge  of 
the  observation  covariance  matrix  is  essential.  Examples  of 
such  applications  include  MVDR  beamforming  and  direc¬ 
tion  of  arrival  estimation  using  MUSIC.  Many  algorithms 
for  estimating  the  covariance  matrix  are  available.  Perhaps 
the  simplest  and  at  the  same  time  most  common  is  given  by 

1  ^ 
m=l 

which  is  the  maximum-likelihood  estimate  given  identical, 
independent,  zero-mean  random  vectors  x^n  with  covari¬ 
ance  R.  Other  estimators  incorporate  information  about  the 
array  geometry.  These  are  commonly  called  structured  co- 
variance  estimators  and  were  introduced  in  [1]. 

When  the  array  changes  shape  significantly  over  an  ob¬ 
servation  interval  the  statistics  of  the  data  vectors  change, 
however.  This  invalidates  the  identical  distribution  assump¬ 
tion  used  to  obtain  (1)  and  the  assumptions  of  most  struc¬ 
tured  covariance  estimation  algorithms.  The  phenomenon 
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of  time-varying  arrays  of  sensors  exists  in  nearly  all  array 
applications.  (No  array  is  truly  time-invariant,  although  they 
may  be  close  enough  to  achieve  the  desired  performance.) 
The  effect  is  exaggerated,  though,  in  towed  sonar  arrays 
which  are  subject  to  underwater  currents  and  the  maneuver¬ 
ing  of  their  parent  platform.  An  array  of  sensors  in  which 
each  sensor  is  mounted  on  a  different  platform  with  its  own 
propulsion  also  constitutes  a  time- varying  array. 

Direction-of-arrival  and  spectrum  estimation  for  time- 
varying  arrays  has  been  addressed  by  a  number  of  authors. 
Direction-of-arrival  estimation  was  addressed  in  [2]  and  [3]. 
Fast  algorithms  for  doing  the  same  which  are  based  on  the 
eigenstructure  of  the  matrix  are  presented  in  [4].  In  [5] 
the  EM  algorithm  is  used  to  estimate  the  power  of  far-field 
sources  using  a  time-varying  array. 

In  this  paper  we  address  the  problem  of  maximimum 
likelihood  (ML)  covariance  estimation  for  time- varying  ar¬ 
rays.  We  proceed  by  defining  a  mathematical  infrastructure 
and  applying  commonly  used  linear  algebra  techniques.  We 
then  propose  several  search  algorithms  to  find  the  covari¬ 
ance  that  maximizes  the  likelihood  under  several  constraints 
imposed  by  the  array  motion.  What  results  may  be  consid¬ 
ered  a  time-varying  structured  covariance  estimation  algo¬ 
rithm. 

2.  DEFINITIONS 

Let  N  be  the  number  of  elements  in  the  array  and  p  n(0  ^ 

be  the  position  of  the  nth  element  at  time  t.  Let  M 
denote  the  number  of  data  vectors  sampled  by  the  array  at 
times  ■**  with  sampling  frequency  Fg.  The 

Tnth  data  vector  we  represent  by  which  is  a  normally 
distributed  complex  random  variable  with  mean  0  and  co- 
variance  R(itn)  =  Rm-  The  time- varying  nature  of  the 
array  implies  that  Rm  need  not  equal  Rm-i-i-  The  problem 
is  therefore  to  estimate  R^  for  m  =  1,  -  ■  *  ,  M. 

We  make  two  assumptions  regarding  the  available  infer- 
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mation.  First,  the  iV  x  1  steering  vector  a(©,  t)  is  known  for 
all  tm  and  for  all  0  €  5^.  The  nth  element  of  the  steering 
vector  is  given  by 


a„(©,f)  =:exp 


A 


(2) 


where  k(©)  is  the  unit  vector  associated  with  the  direction 
0.  Secondly,  the  signal  originating  at  any  direction  is  un- 
correlated  with  signals  originating  at  other  directions.  Also, 
the  sampling  rate  is  such  that  the  sampled  signals  are  in¬ 
dependent  random  variables.  The  time-varying  covariance 
matrix  is  then  given  by 

B^=  [  cr2(e)a(©,t„)a"(e,t„)d©  (3) 

Js^ 

where  (j^(0)d0  is  the  time-invariant  power  of  the  differen¬ 
tial  emitter  at  location  0. 

Since  we  are  interested  in  a  sequence  of  Hermitian  ma¬ 
trices  let  us  introduce  the  following  notation: 


Deflmtion  1  For  positive  integers  N  and  M,  let  Vnm  be  a 
space  such  that  X  6  Vn,m  implies  that  X  =  [Xi ,  •  •  •  ,  Xm] 
where  Xm  €  atUx^  =  X„^. 

Observe  that  we  denote  elements  of  this  space  by  capital, 
bold-faced  letters  with  an  overbar  and  the  mth  element  of 
the  sequence  by  the  same  letter  with  a  subscript.  For  some 
ae  R  and  X  €  we  define  scalar  multiplication  as 

aX  =  [aXi,---,aXM].  (4) 

Similarly,  addition  is  defined  element- wise,that  is  for  X,  Y  e 

Viv,M 

X  +  Y  =  [Xi-hYi,.--  ,Xm+Ym]-  (5) 

It  is  easy  to  see  that  under  these  operations  Vjv,M  is  a  vector 
space  over  R.  For  notational  convenience  we  also  define  the 
following  operations  on  vectors: 

XY=  [XiYi,---,XmYm]  (6) 

X-^=  (7) 

Notice  that  under  vector  addition  as  defined  in  (5)  and  using 
(6)  as  vector  multiplication,  Vn,m  forms  a  non-commutative 
ring  .  The  multiplicative  identity  is  then  the  length-M  se¬ 
quence  of  iV  X  iV  identity  matrices  and  (7)  is  the  multiplica¬ 
tive  inverse  of  X.  With  this  in  mind,  it  would  be  appropriate 
to  refer  to  (6)  and  (7)  as  multiplication  and  inversion  respec¬ 
tively. 

Definition  2  VX,  Y  6  Vn,m  let 

M 

{X,Y)=  ^tr(X«Y„). 

m=l 


We  claim  that  this  is  an  inner  product  on  This  is  a 

result  of  the  following  facts  which  are  easily  proved  for  all 
a  €  R  and  X,  Y,  Z  6 

{X,Y)eR 

(X  +  Y,Z>  =  (X,Z)  +  (Y,Z) 

{aX,Y)  =  a(X,Y> 

{X,Y>  =  (Y,X) 

{X,  X)  >  0  with  equality  iff  X  =  [0,  •  •  •  ,  0] 

Therefore  (V}v,m7  *7 ))  is  an  inner  product  space. 

The  covariance  matrix  sequence  is  an  element,  R,  of 
VnM'  rewrite  (3)  as 

R=  /  a2(0)§(0)d0  (8) 

Js^ 


where 


$(0)  =  {a{0,fi)a^(0,ti),*--  ,a(0,fM)a^(0, W)}* 


The  span  of  $(0)  over  all  0  6  5^  is  a  vector  subspace  of 
VnM'  will  call  this  subspace  Wi,  It  is  clear  from  (8) 
that  R  G  VFi .  Being  a  vector  space,  VTi  is  convex  and  there¬ 
fore  path-connected.  Furthermore,  there  exists  an  orthonor¬ 
mal  basis  for  Wi.  We  will  let  W2  E  Vn,m  be  the  space  of 
all  length-M  sequences  of  non-negative  definite  Hermitian 
matrices.  Since  any  covariance  matrix  is  non-negative  def¬ 
inite,  R  €  IY2.  It  can  be  shown  that  W2  is  also  convex. 
Since  the  desired  sequence  lies  within  both  subspaces  the 
constraint  space  is  their  intersection  W  =  lYi  fl  W2.  As  the 
intersection  of  two  convex  sets,  W  is  also  convex. 

We  remark  that  the  set  of  matrix  sequences  W  may  not 
coincide  exactly  with  the  space  of  matrix  sequences  given 
by  the  model  in  (3),  although  as  has  been  shown  the  latter 
is  a  subset  of  W.  Our  constraint  space  may  contain  elements 
outside  of  the  cone  described  by  (3),  The  discrepancy  be¬ 
tween  the  two  spaces,  and  the  consequences  thereof,  remain 
open  questions. 

We  will  now  derive  the  log-likelihood  function  of  the 
covariance  matrix  sequence  for  the  given  data  set.  The  pdf 
of  the  data  vectors  is  defined  only  for  the  interior  of  1^2 : 

(M 

n 

m-1 

M 

-  Z)  Xm  )  •  (9) 

771=1  / 
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The  log-likelihood  is  then 

M  M 

K'R)  =  -  X)  1^1  ~  X 

m=l  m=l 

M  M 

—  ~  |Rm|  “  (Ry^^Sm)  (10) 

m=l  m=l 

where  we  shall  call 

Sm=Xmx"  (11) 

the  sample  covariance  matrix  at  time  t  m  •  Observe  that  S  ^  = 
and  therefore  the  length-M  sequence  of  all  such  matri¬ 
ces,  S,  is  an  element  of  V}v,m-  We  can  use  the  notation  of 
definition  2  to  simplify  (10): 

M 

l(R;S)=-^ln|iU|-(R-SS>.  (12) 

m=l 

To  find  the  gradient  of  the  log-likelihood  function  we 
will  make  use  of  several  differentiation  theorems  found  in 

Ul 

Theorem  1  For  R,  $  e 

^ln|R-|-x$l  =tr(R-i#). 

Theorem  2  For  R,S,^  e  C^  xjv 

^tr((R+a:#)"^S)  =  -  tr(R-i#R-iS). 

The  directional  derivative  of  the  log-likelihood  along  the 
vector  ^  is 

d  -  -  d  ^ 

— /(R  +  *$;  S)  =  ^  In  |R„  -h 

m=l 

4"  tr((R7Ti  4“  X^yyj)  Sttj.) 

M 

^  y  V  tr(Ry„^$rnRni^Sm) 

m=l 
M 

=  -  53  tr(R-i$„  -  R-iS„R-^$™) 

m=l 

M 

“  ~  {(^rn 

m=l 

=  (R-^SR“^-*R-\#).  (13) 

Therefore  the  gradient  of  the  log-likelihood  function  is  given 
by 

V/(R,  S)  =  R"^SR“^  --  R“X  (14) 


Note  the  similarity  of  this  expression  to  the  analogous  ex¬ 
pression  for  the  gradient  of  the  log-likelihood  in  Burg  et  al 
[1].  Here  the  matrices  have  been  replaced  with  matrix  se¬ 
quences. 

3.  ESTIMATION  ALGORITHMS 

One  possible  estimator  is  the  projection  of  the  sample  co- 
variance  matrix  sequence,  S,  onto  the  constraint  space  W. 
This  is  equivalent  to  selecting  the  point  in  W  that  is  the 
closest  to  S  by  the  standard  distance  metric  for  inner  prod¬ 
uct  spaces: 

d(X,Y)  =  (X-^Y,X~Y)i  (15) 

Because  of  the  similarity  between  this  estimator  and  clas¬ 
sic  filtering  where  a  signal  is  projected  onto  the  subspace 
of  all  signals  which  satisfy  a  certain  constraint,  we  will  re¬ 
fer  to  this  estimator  as  the  sample  sequence  filter.  Since 
W  is  the  intersection  of  two  convex  spaces  we  employ  the 
method  of  projection  onto  convex  sets  (POCS)  in  which  the 
estimate  is  first  projected  onto  Wi  and  then  onto  W^.  This 
iterative  procedure  continues  until  the  improvement  in  like¬ 
lihood  gained  with  each  iteration  is  negligible.  Since  PFi  is 
a  vector  subspace  the  projection  of  a  vector  X  onto  Wi  is 
given  by 

L 

X'  =  53(X,#,)4;  (16) 

l=l 

where  are  the  members  of  an  orthonormal  basis  of  Wi 
and  L  is  the  dimension. 

Projection  onto  W2  for  the  given  inner  product  is  found 
in  [6].  First,  the  eigendecomposition  is  determined: 

x^  =  r^A^r-i  (17) 

Then  the  projection  onto  the  set  of  non-negative  definite  ma¬ 
trices  is  given  by  setting  the  negative  eigenvalues  to  0: 

x;„  =  r.„max(A.„,0)r-i.  (18) 

The  projection  of  the  sequence  X  onto  W2  is  the  element¬ 
wise  projection  of  each  X^^  as  described  by  this  equation. 

The  sample  sequence  filter,  by  its  definition,  finds  the 
sequence  which  is  in  the  constraint  space  and  the  closest  to 
the  sample  sequence  by  the  distance  metric  given  in  (15). 
Experience  has  shown,  however,  that  the  best  estimate  is 
rarely  the  closest  to  the  sample  sequence.  We  therefore  pro¬ 
pose  searching  the  constraint  space  for  the  maximum  likeli¬ 
hood  estimate  using  the  filtered  sample  sequence  as  a  start¬ 
ing  point. 

Each  of  the  search  algorithms  which  we  will  consider 
proceed  by  calculating  a  search  direction,  D  €  V7i,  along 
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which  the  likelihood  function  mus^be  maximized.  That  is, 
in  each  iteration  we  determine  a  D  and  then  find  Ao  such 
that 


Ao 


arg  max/(R 


+  AD). 


Therefore, 

(VUR;  s  -  D),  D)  =  (R-1  (S  -  D)R-1  -  R-^  D) 
=  (VZ(R,S)-R-iDR-\D) 

=  -<R-idR-SD> 


The  updated  estimate  is 


Mi+i)  ^(i) 

R  =  R  +  AoD. 


This  iterative  process  should  be  allowed  to  continue  until 
the  gradient  is  sufficiently  close  to  being  orthogonal  ioWi, 
that  is  until 


(V/(R;S),V/(i;S)) 


Since  D  €  we  know  that  (V/(R;S  —  D),D)  =  0. 
Therefore 


{R-^DR“SD)  =0. 

It  can  be  shown  that  this  implies  that  R“^D^  =  0  for  all 
m.  Therefore  Dm  =  0  for  all  m  which  is  a  contradiction. 

■ 

We  now  concentrate  on  finding  the  direction  which  sat¬ 
isfies  the  necessary  condition  on  the  gradient.  This  is  equiv¬ 
alent  to  finding  D  which  satisfies 


Observe  that  since  I  is  defined  only  on  the  interior  of  and 
D  €  Wi,  the  estimate  will  be  within  the  constraint  space  at 
each  iteration. 

Perhaps  the  most  obvious  approach  to  calculating  D  is 
to  use  the  gradient  in  (14).  D  can  be  the  projection  of  the 
gradient  onto  Wi .  Alternatively,  a  conjugate  gradient  direc¬ 
tion  can  be  calculated  by  incorporating  memory  of  previous 
search  directions.  We  suggest,  however,  a  modification  of 
the  inverse  iteration  algorithm  proposed  by  Burg,  et  al  in 
[1].  Burg’s  algorithm  was  designed  for  estimation  of  a  sin¬ 
gle  matrix  rather  than  a  sequence  of  matrices  but  is  easily 

-  (*) 

generalized  for  sequences.  For  some  estimate  R  select  a 

-(»)  -  - 

search  direction,  D,  such  that  V/(R  ;  S— D)  is  orthogonal 

to  Wi.  Clearly,  if  D  =  0  then  R  is  the  maximum  likeli¬ 
hood  estimate  since  the  likelihood  gradient  is  orthogonal  to 
the  constraint  space  at  that  point. 

Before  the  modified  inverse  iteration  algorithm  may  be 
seriously  considered,  though,  one  must  ask  whether  a  sta¬ 
ble  point  of  the  algorithm  maximizes  the  likelihood  func¬ 
tion  within  the  constraint  space.  That  is,  does  each  iteration 
of  the  algorithm  lead  to  an  improvement  in  likelihood  for  a 
nonzero  search  direction?  The  answer  is  yes,  as  shown  by 
the  following  theorem: 

Theorem  3  Suppose  there  exists  D  ^  0  6  such  that 
Vi(R;  S— D)  w  orthogonal  to  Wi.  Then  there  exists  X£  R, 
A  ^  0,  such  that  /(R  •+■  AD;  S)  >  /(R;  S). 

Proof:  By  way  of  contradiction,  suppose  that 

argmajc/(R-h  AD;S)  =  0. 

This  implies  that  V/(R;  S)  is  orthogonal  to  D.  That  is, 

(V;(R;S),D)  =  0. 


{R-i(S  -  D)R-i  -  R-^  #i>  =  0  (20) 

for  all  L  We  note  that  this  is  a  system  of  equations  which  are 
linear  in  D  and  that  therefore  a  closed-form  solution  exists. 
Since  D  eWi  there  exist  real  Oj  such  that 

i)  =  Y,o‘ih-  (21) 

j 

Substituting  (21)  into  (20)  and  rearranging  we  get 

{R-i(S-D)R-i  -R“^$i) 

=  {R-i(S  -  ^  -  R-\ #i) 

=  (VZ(R;  S),  h)  -  E  a,(R-i'iiR-S  ■ii).  (22) 
j 

Therefore  we  need  only  find  a  6  such  that 

Aa  =  B  (23) 

where 

=  (24) 

Bi  =  (V/(R;S),#i).  (25) 

4.  COMPUTER  SIMULATION 

We  have  simulated  a  uniform  linear  array  (ULA)  consist¬ 
ing  of  A”  =  5  isotropic  sensors  which  is  rotating  with  ro¬ 
tational  velocity  u  about  the  center  element.  The  axis  of 
rotation  is  orthogonal  to  the  axis  of  the  ULA.  There  are  3 
source  fields  impinging  upon  the  array  which  originate  at 
(azimuth,  elevation)  =  (45°,  0°),  (85°,  20°),  and  (110°,  0°). 
Here  azimuth  is  the  the  angle  made  with  the  axis  of  the  ar¬ 
ray  at  <  =  0  and  within  the  plane  of  rotation.  Elevation 
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Figure  1 :  Algorithm  convergence  rate  comparison. 


is  the  angle  made  with  the  plane  of  rotation.  For  example, 
(90°,  10°)  would  describe  a  direction  orthogonal  to  the  ini¬ 
tial  array  axis  and  10°  above  the  plane  of  rotation.  Each 
of  the  sources  are  assumed  to  be  narrowband  with  wave¬ 
length  A  and  the  separation  between  elements  in  the  ULA 
is  The  power  of  each  source  at  the  array  is  30dB,  15dB, 
and  20dB  respectively.  Receiver  noise  is  OdB.  The  rota¬ 
tional  velocity  of  the  array  is  a;  =  27r  rad/sec,  the  sampling 
rate  is  =  32s"  ^  and  the  number  of  samples  collected 
is  M  =  IG.  Therefore,  the  array  gathers  16  data  vectors 
while  completing  a  half  rotation.  Since  the  statistics  of  the 
data  vector  change  dramatically  over  this  observation  inter¬ 
val  one  expect.'*  the  covariance  estimator  in  (1)  to  perform 
badly.  It  is  unclear  even  what  steering  vector  to  use  with 
this  covanance  estimate  in.  for  example,  an  MVDR  or  MU¬ 
SIC  estimator  7>ia!  makes  this  scenario  a  good  candidate 
for  covanance  matnx  sequence  estimation. 

Each  of  the  algonthms  considered  begins  by  calculating 
the  sample  co\anancc  matnx  sequence  and  filtering  it  by 
the  method  of  POCS  using  the  projections  in  (16)  and  (18). 
The  arc  obtained  b>  Gram- Schmidt  orthonormalization 
of  the  set  of  vectors  4^(9)  w  here  0  is  a  discretization  of  the 
2-sphere.  For  this  example,  the  dimension  of  Wy  is  L  = 
76.  Note  that  the  dimension  of  is  N*^M  =  400  and 
we  have  managed  to  eliminate  324  parameters  of  the  matrix 
sequence  by  applying  information  about  the  motion  of  the 
array. 

The  ML  search  routine  follows  the  filtering.  We  applied 
a  gradient  and  conjugate  gradient  search  algorithm  in  addi¬ 
tion  to  the  modified  inverse  iteration  algorithm.  The  likeli¬ 
hood  of  the  estimate  at  each  iteration  is  plotted  in  Figure  1 
for  each  of  the  algorithms. 

Observe  that  the  gradient  and  conjugate  gradient  algo¬ 
rithms  converge  to  points  with  the  same  likelihood.  The 


Figure  2:  MVDR  spectrum  estimated  from  the  first  matrix 
in  the  sequence.  The  dashed  line  is  that  obtained  from  the 
sample  sequence  filtering  procedure.  The  solid  line  was  cal¬ 
culated  from  the  ML  sequence. 


conjugate  gradient  reaches  this  point  in  fewer  iterations,  which 
is  to  be  expected.  However,  the  convergence  point  of  the  in¬ 
verse  iteration  algorithms  exceeds  the  likelihood  of  the  es¬ 
timate  obtained  from  either  gradient  algorithm  after  only  a 
few  iterations.  Inspection  of  the  likelihood  gradient  at  what 
appears  to  be  the  convergence  point  of  the  gradient  algo¬ 
rithms  reveals  that  it  is  not  orthogonal  to  Wi  and  that  while 
successive  iterations  yield  only  slight  improvement  in  like¬ 
lihood,  they  have  failed  to  reach  a  local  maxima.  One  possi¬ 
bility  is  that  they  have  stumbled  upon  a  ’'ridge’*  in  the  like¬ 
lihood  function.  It  is  clear  that,  in  this  example  at  least,  the 
inverse  iteration  algorithm  reaches  a  solution  in  fewer  itera¬ 
tions  than  even  the  conjugate  gradient  algorithm.  It  should 
be  noted,  however,  that  finding  the  solution  to  (23)  requires 
more  compuation  than  calculating  the  likelihood  gradient 
and  projecting  it  onto  Wi . 

To  demonstrate  the  validity  of  the  ML  estimate,  the  MVDR 
spectrum  corresponding  to  the  first  matrix  in  the  sequence 
has  been  calculated  and  plotted  in  Figure  2,  The  spectrum 
is  calculated  using 


_ 1 _ 


(26) 


Ri  is  the  first  matrix  in  the  sequence  obtained  by  the  in¬ 
verse  iteration  algorithm  since  the  other  two  algorithms  failed 
to  produce  a  maximum  likelihood  estimate.  The  position  of 
each  of  the  sources  is  easily  ascertained  from  the  plot  as  is 
a  feeling  for  their  intensities.  Also  plotted  is  the  spectrum 
obtained  from  just  the  sample  sequence  filter.  While  peaks 
which  correspond  to  two  of  the  sources  can  be  seen,  the 
third  is  lost  and  the  background  noise  is  quite  high.  This 
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demonstrates  the  necessity  of  the  ML  search  algorithms. 


5.  CONCLUSIONS 

We  have  developed  an  algorithm  for  estimating  the  sequence 
of  matrices  which  are  the  covariances  of  the  data  vectors  of 
a  time-varying  sensor  array.  Since  each  matrix  in  the  se¬ 
quence  is  strucmred  to  satisfy  (3)  this  may  be  classified  as 
a  structured  covariance  estimation  algorithm  in  which  the 
sequence  of  matrices  itself  is  structured.  This  method  will 
have  good  performance  for  time- varying  arrays  in  which  the 
array  motion  is  periodic,  as  with  the  rotating  ULA,  since 
the  constraint  space  basis  need  not  be  continuously  recalcu¬ 
lated,  which  relieves  the  processor  of  some  of  the  compu¬ 
tational  burden.  It  has  been  demonstrated  that  the  modified 
inverse  iterations  algorithm  can  converge  faster  and  more 
reliably  than  a  simple  gradient  search  algorithm,  although 
with  increased  computational  complexity  per  iteration. 
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